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Abstract 

We propose a quite general model of active media by consideration of the inter- 
action between pacemakers via their phase response curves. This model describes a 
network of pulse oscillators coupled by their response to the internal depolarization 
of mutual stimulations. 

First, a macroscopic level corresponding to an arbitrary large number of oscilla- 
tory elements coupled globally is considered. As a specific and important case of the 
proposed model, the bidirectional interaction of two cardiac nodes is described. This 
case is generalized by means of an additional pacemaker, which can be expounded 
as an external stimulater. The behavior of such a system is analyzed. Second, the 
microscopic level corresponding to the representation of cardiac nodes by one- and 
two-dimensional lattices of pulse oscillators coupled via the nearest neighbors is de- 
scribed. The model is a universal one in the sense that on its basis one can easily 
construct discrete distributed media of active elements, which interact via phase 
response curves. 
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1 Introduction 



Representation of an active distributed system by ensembles of coupled excitable or oscil- 
latory elements is very useful method of the analysis because it allows to understand main 
dynamical processes inherent in the considered medium. As is known, this approach goes 
back to the model of Wiener and Rosenblueth [Wiener & Rosenblueth, 1946], according to 
which a medium consists of single elements being in one of three possible states: excited, 
refractory or rest. Later such models as coupled limit cycle oscillators and chaotic maps 
[Kaneko, 1990; Shibata & Kaneko, 1998] have played an important role not only in a quite 
realistic description of active media but also in the understanding of a possible behavior of 
systems far from equilibrium. Many useful concepts like phase-locked patterns, synchro- 
nization and spatio-temporal chaos have become popular due to detailed studies of similar 
nonlinear models [Kuramoto, 1984; Kuramoto, 1995; Winfree, 2000]. 

Investigations of such an example of an active medium as cardiac tissue are of signif- 
icant scientific interest owing to vital importance of its rhythm stability. Real heart cells 
exhibit oscillatory properties (can be reset and entrained), they arc excitable and have 
a refractory time, during which they do not respond to external stimulation. Hence, the 
heart can be considered as consisting of oscillatory (conductive cardiomyocytes, which have 
automaticity) and excitable (contractile heart cells, which do not initiate electrical activity 
under normal conditions) elements. 

Due to extraordinary complexity of the heart, many qualitative discrete and continuous 
models have been tested. Majority of computational models of cardiac tissue of last gener- 
ation takes into consideration the kinetics of excitable cells, how the excitation propagates 
from cell to cell and how contractile cardiomyocytes are arranged and connected in space 
[Clayton, 2001]. Such models mainly serve for studying sustained by re-entrant activity 
lethal arrhythmia - ventricular fibrillation, during which the spatio-temporal behavior is 
very complex [Clayton et al., 2006]. 

Other models treat the cardiac tissue as an active conductive system, taking into ac- 
count oscillatory properties of heart cells. In this case the cardiac rhythms can be described 
on the basis of the dynamical system theory (see e.g. [Courtemanche et ai, 1989; Gold- 
berger, 1990; Bub & Glass, 1994; Glass et a/., 2002; Loskutov et aL, 2004] and refs. therein). 
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Hereinafter we hold this approach. 

Under normal conditions the electrical activity of the heart (action potentials) is spon- 
taneously initiated in a region of the right atrium, sino-atrial (SA) node, so-called lead- 
ing pacemaker. Automatic excitation is a distinctive feature not only of the cells of the 
SA node, but also of other conductive heart cells, so-called latent pacemakers. In addi- 
tion, contractile cardiomyocytes can initiate a spontaneous action potential in pathology. 
Electrophysiological studies have suggested that the activity of cardiomyocytes with auto- 
maticity (e.g. P-cells of the SA node, of the atrioventricular (AV) junction, Purkinje cells) 
can be modulated by current pulses stimulating (super-threshold depolarizing) applied ex- 
tracellularly [Jalife & Moe, 1976; Sano et al., 1978; Jalife et ai, 1980; Antzelevitch et ai, 
1982]. 

Effects of external stimuli on biological oscillators are observed in a wide range of 
species. Experimentally obtained characteristics can be represented by a phase response 
curve (PRC) [Jalife & Moe, 1976; Antzelevitch et al, 1982; Reiner & Antzelevitch, 1985]. 
To establish the shape of the PRC experimentally, stimulation of an oscillator at various 
phases of its intrinsic cycle is applied [Jahfe et al, 1980; Guevara & Shrier, 1987]. It 
has been found that in different pacemaker cells early stimuli delay the next pacemaker 
discharge and late pulses advance it. Therefore, the typical PRC shape is bipliasic [Jalife 
et al, 1980; Reiner & Antzelevitch, 1985]. 

The rhythm of autonomous biological oscillators can also undergo an external periodic 
perturbation (e.g. activity of cells of the AV junction is subjected to sinus rhythm), 
depending on both the stimulus magnitude and its phase within the cycle. It is known that 
when the frequency and the amplitude of the external periodic stimulation are varied, a 
diversity of phase diagrams can be established between the stimulus and the self-sustained 
oscillator (see e.g. [Loskutov et al, 2004]). In some situations the rhythm of the biological 
oscillator is entrained (or phase-locked) to the external stimulation so that for each M 
cycles of the stimulation there are cycles of the autonomous oscillator rhythm. This 
occurs at a fixed phase (or phases) of the stimulus and is called M : N phase-locking or 
entrainment, which appears as a time-periodic sequence. In particular, entrainment of 
1:1, during which the rhythms of the oscillator and external stimulus are matched, is 
defined as synchronization phenomenon. 
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In the present paper we develop a general simplified model describing a network of pulse 
oscillators coupled by their response to the internal depolarization of mutual stimulations. 
Our primary aim is to keep the model as simple as possible and to introduce a minimal 
number of parameters. Therefore, in our model the pacemakers are fully characterized by 
their intrinsic cycle length and are represented as pulse oscillators. Their interaction is 
described by PRCs. At first, we will consider two bidirectionally interacting pacemakers to 
demonstrate the basic concepts of the model. Then we will apply this approach to construct 
a pacemaker network model with global coupling. As the following step, we will analyze 
two specific cases of this PRC based model of coupled pulse oscillators: two and then three 
interacting cardiac nodes. An additional pacemaker can also be expounded as an external 
stimulater. Our further intention is to go on to the next (microscopic) level and represent 
each pacemaker as an ensemble of interacting oscillatory elements. Extrapolation of our 
approach to one- and two-dimensional matrixes (or lattices) of pacemaker cells interacting 
via nearest neighbors concludes the present study. 

2 Development of the General Model 

In this part we construct a system with two pacemakers and then consider a quite general 
model of a set of interacting pacemakers coupled by their PRCs. 

2.1 Two interacting pacemakers: outline of the approach 

Consider two interacting pulse oscillators (or pacemakers) A and B with intrinsic periods 
of their autonomous beating Ta and respectively. An interaction between oscillators is 
governed by so-called phase response curve (PRC). This means that a phase shift of one 
of the oscillators happens as a result of an impact of the another one. To construct an 
adequate mathematical model, it is necessary to accept some restrictions concerning the 
character of the interaction. We describe them briefly [Ikeda, 1982; Glass et al, 1986; 
Glass & Zeng, 1990]. 

1. The phase of the disturbed oscillator is shifted to a new value instantly after an 
impact. 

2. The phase shift depends only on two main parameters: a) on the phase difference 
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of oscillators and b) on the influence strength. In turn, this influence strength depends on 
its amplitude and the coupling coefficient of the oscillators. In a real system the coupling 
coefficient is the average factor that shows how the strength of the pulse decreases during 
its passing from one oscillator to another. Thus, the phase shift $ determining a new phase 
of the disturbed oscillator with period T can be represented as follows: 

$ = A/T = $(<^,£), (1) 

where A is the time shift of the disturbed oscillator, ip is the phase difference of the 
oscillators and e is the influence strength. 

Pacemakers can be represented as a set of separated firing peaks on a time scale. Assume 
that the instants of the last firings of the oscillators A and B are a and h respectively 
(Fig. Note that a and h are the moments of the impacts after all previous phase shifts 
of the oscillators. In other words, one can observe the oscillators' firings at these particular 
instants. Then it is necessary to analyze two cases: 

1. b < a. This is the case 1 in Fig. ^ i.e. when the oscillator B has fired before the 
oscillator A. Let us follow the dynamics of the system in real time. The nearest 
event, that affects on the further behavior of the entire system, is an appearance of 
the pulse A at time a. Let us stop on at this moment and make the forecast. To this 
end we define a concept of the moments of expected firings of the oscillators, i.e. a'^ 
and b'^. Let us imagine that we have shifted back in time with respect to the moment 
a. Since the oscillator A has not fired yet, one should expect the appearance of the 
next A and B pulses at the moments = a and b'^ = b + Tf, respectively, where Ti, 
is the period of B. We call this situation as "A fires and B is at an expected state" 
and denote symbolically as {a,b^). Now we consider the instant a. Since A fires, the 
next expected values can be transformed to 

Kext = b + Tb + Ab{(fib, eb) = b^ + Ab{ipb,€b), 
where Ab{ipb,£b) is the time shift of the oscillator B due to the impact of A. It 
depends on the phase ipb of the pacemaker A with respect to B and the influence 
strength Eb- The phase ipb can be calculated as follows 



or, in terms of the expected values, 



a' 




(mod 1). 



Note that the phase ipb is a positive value, and it belongs to the [0, 1] segment (neg- 
ative values in the two previous expressions are eliminated by the mod 1 operation). 

To determine which oscillator fires next, one should compare a^g^^ and &^ext- If 
^next < Kexty then A fircs and B remains at an expected state until b^^^^, i.e. the 
system moves to the state (a, V)next- Otherwise, if b'^next < 'Kiexv then B fires, and A 
jumps to an expected state and the entire system's state becomes (a*^, b)next- 

2. b > a. This is the case 2 in Fig.^ This inverse situation is analogous to the previous 
one with the difference in speculations owing to the firing of the pacemaker A prior 
to B. Then expected values a*^ and b^ can be written accordingly: = a + and 
b^ = b. One can call this "i? fires and A is at an expected state" and denote 

by (a^, b). The next expected values are given by the following expression: 



where Aa{(pa,£a) is the time shift of the oscillator A. It depends on the phase of 
the oscillator B with respect to A, i.e. ipa = (^^ — 0'^)/Ta (mod 1), and the influence 
strength Sa- Further analysis is also similar to the case 1. Namely, if fo^g^.^ < a^ext; 
then B fires and A jumps to the expected state a^g^,^, i.e. the system moves to 
the state {a^,b)next- Otherwise, if a^ext < Kext^ then A fires and B remains in the 
expected state and the system state becomes (a, b^)next- 

Summarizing the above calculations, the model can be represented by the following 



In the notions of the expected values, which we denote for convenience as = aW = b, 



next 



'next 



a" + Aa{(Pa,ea), 

b' + n, 



scheme: 




6 



the dynamics can be described by the following difference equation: 



bn+1 J \bn 



, ttn <hni and then A fires at time a„ 

(3) 

(v^n) ^a) | ^ '^^ < ^^(1 then B fires at time 



where: 



\ " II "'■II j , , !, ( I 

= ^ ^ ^ (mod 1) , = ^ ^ ^ (mod 1) 



To simulate the dynamics of two interacting oscillators coupled by PRCs, it is necessary 
to carry out the iterative process Q for the expected pulses and put sequentially in the 
time scale the firing moments of the oscillators A and B depending on the result of the 
comparison Un and 6„. An investigation of the case of two interacting pacemakers in more 
detail with the description of the possible modes of behavior is given in Section 3.1. 

2.2 Derivation of the basic model equation 

Assume that there are autonomous pulse oscillators (or pacemakers). Suppose that 
all the pacemakers are different. This means that each has its own intrinsic cycle length 
Tj, i = 1, . . . ,N, and the beating amplitude. To define coupling between pacemakers, it 
is necessary to determine the topology of the system space. In other words, one should 
specify the nearest neighbors of each pacemaker in the space. Vice versa, it is obvious that 
the determination of coupling between the elements of such spatially discrete system sets 
its topology. 

First of all we develop the general model of A^ mutually coupled pulse oscillators. 
Suppose that all the pacemakers interact with each other, i.e. so-called global coupling 
is realized. The model derived in Section 2.1 can be easily generalized to the case of A^ 
pacemakers. We operate with the expected values introduced in Section 2.1, the real firings 
of the pacemakers are found by the analysis of the expected impacts series. 

Let the set of expected firings {aj}i^...^Ar be located in a time axis (see Fig. This 
means that in the absence of coupling, pacemakers strike at these instants. Suppose now 
that some oscillator acts on another one by means of the PRC {Aij{(f^\eij)}i^,,,^N, where 
is the phase of the j-th pacemaker with respect to the i-th one and Sij IS db total 
parameter defining coupling between the j-th and i-th elements. 
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The next values of the expected firings can be calculated by the same manner as in 
Section 2.1. Because the j-th oscillator appears before all others (see. Fig. 12)), it does not 
undergo any influence and fires as a real impact of the j-th pacemaker. Thus, the given 
event makes the shifts of all other oscillators according to the set of the PRCs {Aj^}. The 
j-th oscillator is shifted to a new expected moment as an unperturbed one, i.e. by adding 
its own cycle length Tj. To get the next sequential expected values, one should make the 
same procedure with the newly obtained expected pulses. The dynamics of the system can 
be easily represented as the following iterative relation: 

f T- i = j 

K+i = «n + I ^^ij^ ^ ^ _^ J ■ < = min«}i=i,...,7v, (4) 



where 



= \^ (mod 1) 



It is convenient to rewrite the PRCs by normalizing {Ajj} on the intrinsic pacemaker 
cycle lengths Tj and to define as follows: 



A,, = ■ T,}, e [0, 1] 



where {fij{^pli,£ij)} are the dimensionless functions, which are also called the PRCs. The 
real pacemakers can have identical nature but differ in the intrinsic cycle lengths. Then 
the form of dimensionless PRCs is identical for them, i.e. their f{(f,e) coincide, while 
A{(f,e) are different. Moreover, it is convenient to use the functions {fij{fli,£ij)} in the 
construction of the equations for the dimensionless phase differences between pacemaker 
pairs. In Section 3 we will demonstrate this approach for the systems of two and three 
coupled pulse oscillators. 

3 Applications 

3.1 Two interacting pulse oscillators 

Let us investigate the system of two interacting pacemakers in detail. As well as in Section 
2.1, suppose that the system consists of two oscillators A and B coupled by means of the 
phase response curves Aa((^°,ea) and A;,((^|^, It is clear that the map Q specifying 
the dynamics of such a system is a particular case of the general model (0)) restricted to 



N = 2. Rewrite Q using the expressions for the dimensionless PRCs: 

cin < bn, and then A fires at time a„ 



bn+1 j \bn 



Ta 



+ < 



(5) 



b 



bn <an, and then B fires at time bn, 



where: 



-La J-b 



Let us assume that the pacemakers have one the same nature. Hence, one can accept 
fa{ip,e) = fb{'^,e) = f{ip,e), where f{^,e) is an one-parametric function. The parameter 
e integrally determines a total influence of one oscillator on another one. In the symmetric 
case we have ea = Sb = s. In our further consideration we do not accept this symmetry. 

From the mathematical point of view the system is a map of the real plane into 
itself depending on four parameters: T^, T^, Ea and e^. It is easy to comprehend that the 
effective parameter that changes the behavior of the system (jSj) is the dimensionless ratio 
5 = Th/Ta- Specifying parameters Ea, Eb, S and the function f{(f, e), one can iterate directly 
the equations (jH)) marking on the time axis the values of the real firings of the pacemakers 
A and B. Some examples of this simulation are presented below. 

In spite of the fact that the expression defines the two-dimensional map, the values 
an and bn increase gradually because they represent the time series of the expected firings 
of the oscillators. Therefore, trajectories of the map (0) are infinite. Thus, this is not 
informative for us. It is more important to derive a map that determines the dynamics of 
the phase difference of the pacemakers. 

Let us introduce the dimensionless phase difference of the pacemakers A and B: 

(^n bn 

a 

The choice of Ta as a normalization factor is not essential. Having selected Tf,, as a result 
we obtain the similar expressions. Subtracting the second equation of the system Q from 
the first one and dividing the result by Ta we get the following: 

Xn^, = {^- + ^-^f{^ (^Od 1), , Xn<0, 

Xn + f (Xn (modi), Ea) - 6, X„ > 0, 
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where 6 = T^/Ta- Here we take into account that < if < and x„ > if a„ > 6^. 

Let us make a brief analysis of the developed map. Because in general x G {~oo; oo), the 
equation (jHl) represents a one- dimensional nonlinear map of the real axis into itself. Note 
that the map can not be reduced to a circle map by the restriction of x to the range [0; 1] 
as it is usually done for two pacemakers interacting by PRCs. It is essentially asymmetric 
with respect to changing x to —x (see Fig. IS)). If f{x,e) is a nonmonotonic function, then 
the map © is nonlinear, and it can exhibit a big variety of the behavior: from complex 
periodic motion to chaotic dynamics. Owing to x G {—oo; oo), in a rigorous sense it is not 
a difference of phases of the pacemakers. In this context one can call x as the generalized 
phase difference. Analyzing Eq. © one can find out which oscillator, A or B, fires at the 
given discrete time n. This depends on the sign of x: A fires if x„ < and B fires when 
Xn > 0. Thus, it makes possible to determine the phase- locking degree of the pacemakers. 
However, taking into consideration only the values x„, we may not reconstruct the initial 
time series of firing events of the pacemakers A and B. Using Xn, one can say only about 
their phase difference. 

Let us show how the model equations (0) can be applied to the investigation of the 
behavior of two interacting pacemakers. 

Analysis of real systems [Glass et ai, 1987] shows that the function f{x,e) can take 
different forms. But, as a rule, it obeys a number of general properties. For example, 
f{0,e) = f{l,e) = 0. Usually it has a maximum and a minimum. Sometimes instead of 
extrema it has breaks. Let the function f{x, e) be in an elementary form that is often used 
(see, e.g. [Glass & Zeng, 1990]). Namely, let f{x,e) = esin27ra;. This leads to an array of 
dynamical equations: 



K+l I \bn 



+ < 



e„sin(27r(^»)T, 



bn < an- 



J x„ + 1 - fe;,sin(— ^ (modi) ■ 27r), x„ < 0, . . 

Xn+l — \ [O) 

Xn + Ba sin(x„ (mod 1) • 27r) —6, Xn > 0. 
In Fig. 01 examples of direct simulation of the system ((7j), (jH)) are presented. In the 

left column some possible phase-lockings found on the basis of Eq. ^ are displayed. In 
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the right column the corresponding map (jH)), its periodic orbit and values of Lyapunov 
exponents are shown. Fig. Ell represents the existence of the chaotic behavior for the 
system at Sa = 0.1205, Eb = 1.2845 and 6 = 0.6465. 

As a conclusion of this section it should be noted that the system of two bidirectionally 
acting pacemakers has been already intensively investigated (see, e.g., [Ikeda et ai, 2004] 
and references herein) as a qualitative model of cardiac arrhythmia known as modulated 
parasystole. The interacting oscillators represented the pairs of cardiac pacemakers such 
as: the sinoatrial (SA) node and the ventricle contracting by different factors, e.g. ectopic 
pacemaker, etc. Hereby the authors used various kinds of PRCs f{x,e) approximating 
the experimental data on stimulating the cardiac cells of animals by single electric current 
pulses. Recently in the paper [Loskutov et ai, 2004] the system of two interacting pacemak- 
ers similar to ((Zj), (jH)) has been analyzed. In this paper taking into account the refractory 
time the various types of the smooth PRCs were examined and the bifurcation diagrams 
of possible phase- lockings were constructed. However, in [Loskutov et ai, 2004] the be- 
havior regimes when the firings of the pacemakers are not alternated, i.e. cases of strong 
discrepancies in the intrinsic cycle length (T^ ^ T5 and vice versa), were not investigated. 
The system (j?)), (jH)) is more general and takes into account all possible variants. 

3.2 Three pacemakers 

Below we explain why the special case of three interacting oscillators is worth individual 
attention. First, as is known there are three drivers of the rhythm in the cardiac conductive 
system: the SA node is a leading pacemaker, the AV junction and Purkinje fibers are 
latent pacemakers, which under normal conditions are suppressed by the sinus rhythm. 
However, at violations of conductivity and pulse initiation, the cardiac pacemakers can 
influence on each other, i.e. so-called bidirectional coupling may be realized. Second, in 
pathology a group of contractile cardiomyocytes can also initiate action potentials: an 
ectopic (abnormal) cardiac pacemaker may emerge and start to compete with the SA node 
for leading the heart rhythm. Third, at stimulating the cardiac tissue by external current 
impulses (cardio-stimulation), the cardiac rate is changed. The external stimulaters can 
be naturally included in our general model of interacting pacemakers as additional leading 
pulse oscillators. 
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Thus, representation of the heart conductive system as at least three coupled au- 
tonomous oscillators (Fig. 0} is very useful for understanding which influence of an ad- 
ditional pacemaker exerts on a system of two bidirectionally interacting drivers of the 
rhythm (i.e. the SA and AV nodes) considered above. Investigation of the possible behav- 
ior modes of a larger amount of interacting pacemakers turns out to be a very complicated 
both analytical (mathematical) and numerical problem. For example for a system of five 
coupled pulse oscillators we have 25 various functions and 29 different parameters (25 val- 
ues of Eij and 4 independent ratios Ti/Tj). This is due to the fact that, in general, the 
cardiac pacemakers have various frequencies (or intrinsic cycle lengths Tj ) and the different 
nature. This means that the PRCs {fij{x,eij)} determining coupling between a pair of 
pacemakers have different forms. For some heart pacemakers PRCs have been measured 
by the direct experiments on cardiomyocytes of animals (see [Glass et ai, 1986]). Other 
PRCs can be chosen using general principles based on the nature of nodes or on the basis 
of the collateral measurements [Ikeda et ai, 1988]. 

Applying the general model equations (jlj) for three interacting pacemakers, one can get 
the following system: 

/T. \ 







f an \ 


bn+1 




bn 









if a„ < 6„, c„; A fires at time a„ 



V ^cb {Vn^^cb) J 



if bn < (in, Cn, B fircs at time 6„ 



(9) 



Abe [^i',e,,) 



iicn<an,bn] C fires at time c„ 



where 



ba 



ab 



an 


- bn 






bn 


- an 




T 


Cn 


an 


T 



(mod 1) , 
(mod 1) , 
(mod 1) , 



(Xn Cn 



n 



(mod 11 



(mod 11 



fmod 11 



The response functions Aij are supposed to have a form Ajj((y9^, Eij) = fiji^pl^, Sij)Ti, y?^ G 
[0; 1], i,j = a,b,c. The example of a system of three bidirectionally interacting cardiac 



12 



pacemakers: the SA node, the AV junction and ectopic pacemaker, is shown in Fig. 

It is convenient to study Eqs. Q introducing phase differences of the pacemakers. Let 
us define, 



cifi bji 



Xr, 



Vn 



Cfi bj-i 



T 

-J- n 



a 



Tr 



Then for the phase differences we obtain the map of the real plane into itself: 



Vn+l 



Vn 



+ 



f a- fba {{Xn},eba 
I I ■^n Vn 



+ < 




(10) 



Xn> 0, yn> 



Vn < Xn, yn<0 

P - he {{yn] ,£bc 

Here we denote {x} as the fractional part of x. The conditions in the right hand side 
of (Uni) divide the pi ane into three areas corresponding to the firing pacemaker (see Fig. EI). 
The map ()1()|1 has breaks on the boundaries of these areas. 

Suppose that the pacemaker B is an external stimulus (see Fig. |3Jd). Then the ex- 
pressions dn}, pO|) become simpler. The stimulus B acts on the pacemakers A and C but 
does not experience any influence. In this situation fba{x,eba) = fbc{x,ebc) = and the 
model (P), (fTUj) can be rewritten as follows: 



r f Ta 





/ 






{ an \ 




bn+1 




bn 


V 


Cn+1 / 







+ < 



( A,,( 
Tb 



^n.eab] \ 



li Un < bnyCn] A fircs at time a„ 



if hn < Cn] B fircs at time 6„ (11) 



Aae«^£ac) 

I , if c„ < a„, bn] C fires at time Cn 

Tr 
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The corresponding expression for the phase differences takes the form: 



Vn+l 



+ 



/ a 



+ < 




(12) 



Xn>0, yn>0 



Un < Xn, Vn < 0. 



In Fig.Elsome examples of phase- lockings and complex dynamics for the system 0, (fTn|l 
are shown. The PRCs are chosen in the sinus form, i.e. fij{x, Sij) = Sij sin(27rx). Tj and Sij 
are indicated in the legends. It is obvious that the model of three interacting pacemakers 
requires further investigation in detail. In particular, we believe that appropriate choice 
of influence strength and period of external stimulation in the equations ()12|) will remove 
the system from undesirable complex behavior to the normal heart rhythm. It will be 
represented within forthcoming works. 

4 Lattices of Coupled Pulse Oscillators 

Let us demonstrate briefly the way how one can approximate discrete distributed media on 
the basis of the general model of coupled oscillators Q. Considering a cardiac pacemaker 
on the microscopic level, one can represent it as a large group of cells, which generate the 
heart rhythm and synchronize their action potentials to initiate the cardiac contraction. 
Thus, instead of examining a single pacemaker we should construct a lattice of coupled 
pulse oscillators. In this paper we restrict ourselves to the one- (chain) and two-dimensional 
(lattice) cases. 

Let us suppose that the autonomous pacemakers are located in sites of the two-dimensional 
square lattice of the size N x M. An element of the lattice with coordinates (i,j) is des- 
ignated as A*'-', where i = 1, . . . ,N and j = 1, . . . , M. We restrict consideration to the 
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homogeneous medium and accept a number of limitations on anisotropy. This means that 
the pacemakers of the lattice are identical, i.e. they have the same cycle length Tjj = T, 
i = 1, . . . , N, j = 1, . . . , M (however, in fact, cells from the periphery of the sinus pace- 
maker show the shortest cycle length, although the centre acts as the leading pacemaker 
site) . This limitation decreases a quantity of parameters of the system and hence simplifies 
the study of the model. Now we should define coupling between elements. In works devoted 
to the coupled map lattices (CML) two main types of coupling are usually assumed: the 
nearest neighbors and global couplings (see, e.g., [Kaneko, 1989]). Since in the previous 
sections we supposed that pacemakers interact with each other, this time as an example we 
consider lattices with the coupling to nearest neighbors: first, a two-dimensional lattice, 
and then a chain of coupled pulse oscillators. 

Within the square lattice each pacemaker A^'^ interacts with four of the neighbors 
according to the schematic picture in Fig. IHK- Taking into account the limitation of homo- 
geneous medium, one can assert that all couplings of the lattice are identical, i.e. each two 
adjacent elements interact with each other by a general law defined by the identical PRCs 
/(x, e). Moreover, suppose that the coupling between a pair of elements is isotropic in such 
a sense that = equal to one of the values ei or 62 depending on 

the relative position of elements. This means that there is an anisotropy of the influence 
strength in vertical and horizontal directions. In other words, if two pacemakers are the 
neighbors in the vertical direction, they interact by f{x,ei), and if they are horizontal 
neighbors, then they are coupled by /(x, £2)- 

Note that all these restrictions are made for the simplification of the analytical form of 
the resulting model. In general, it is possible to write the expressions for a two-dimensional 
lattice of coupled pulse oscillators without any limitations. It is a subject of a standalone 
investigation and lies beyond the framework of the given study. 

Let us write equations that determine the iterative dynamics of the expected firings of 
the pacemakers \a'^'Hi=i.....N on the basis of the approach represented in Sec. 2. 2. To obtain 

j = l,...,M 

the (n+ l)-th value of the individual element a*'-' , it is necessary to analyze all elements of 
the lattice since they are coupled with each other by means of a local coupling. In other 
words, the considered element cannot be affected by others at the n-th time step because 
the latter is suppressed by the influence of other elements and remains at the expected state 
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until the {n + l)-th step. Then dynamics of the model can be described by the following 
expression: 



-I, J 



/(<^(..)(o,+i), 
/(^(..)(o,-^ 
/(^(..)(+i.o), 
/(^(..)(-i,o), 

0, 
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gmin 
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(Jj 
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if a-i'^- = 
otherwise 
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where the phases are 
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The constructed model demands detailed investigation on the basis of the approach 
developed for the CML (see, e.g., [Kaneko, 1989]). It will be represented in the succeeding 
works. 

As the second example let us consider a chain of the identical pulse oscillators coupled by 
the nearest neighbor principle. We restrict ourselves to a homogeneous case with anisotropy 
of the right and left direction in the influence strength between the nearest neighbors. The 
schematic picture of the chain is shown in Fig. ^p. Similarly to the above consideration 
one can get: 

if at = a™''', 



''n+l 



T, 



f{^li+\e2)-T, if a: 

f{^ll-\e,)-T, ifajr^ 



i+l 
n 



min{a^}i 



=l...Ar 



(14) 



0, 



where the phases are the following 



T 



otherwise 



T 



If ei = (or £2 = 0), then Eqs (fT^ define the so-called open- flow model [Willeboordse & 
Kaneko, 1994]. 
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Because in this work we present a general approach of developing models without the 
detailed analysis of their behavior, the type of boundary conditions for both lattices has 
not been indicated. Hence, to investigate such systems analytically or numerically, one 
should set the boundary conditions along with the PRCs f{x,e). Usually the boundary 
conditions are chosen as periodic, i.e. a^'^^'^^ = a*'-'; W^^'^ = a^'^ for the two-dimensional 
lattice and 2*"*"^ = a* for the one-dimensional one. For the open-flow model a condition of 
the fixed left boundary, = const, is frequently accepted. 

The described models (0}, ()13|) and (fT^ admit generalization to a natural inhomoge- 
neous case by placing different intrinsic cycle lengths of the pacemakers, PRCs and influ- 
ence strengths for various groups of elements. However, consideration of inhomogeneous 
anisotropic lattices is extremely difficult problem even for numerical analysis. The first 
attempts of investigating inhomogeneous lattices of coupled maps (ICML) were described 
in [Vasil'ev et al, 2000; Loskutov et a/., 2002; Rybalko & Loskutov, 2004]. 

5 Summary and Limitations 

In the present study we propose a quite general discrete model of active media by in- 
troducing a simple phase response curve interaction between leading centers. We have 
shown that the PRC can be a useful "tool" for representation of the interaction between 
pacemakers in cardiac tissue both on a large and small scales. This PRC based model to- 
gether with demonstrating complex (chaotic) behavior, can describe the entrainment and 
synchronization phenomena of interacting pulse oscillators. It can also aid to understand 
their response to an external stimulus with variable intensity and duration (see Fig. HJd), 
as previously observed in experimental studies [Jalife et ai, 1976; Jalife et ai, 1980]. 

Starting with consideration of two interacting pulse oscillators and introducing new 
concepts of expected values, we have extrapolated our PRC based approach to investigate 
the mutual influence among an arbitrary large ensemble of pacemakers. The specific cases 
of the proposed model show that it can be very useful for investigating the dynamical 
interaction of cardiac nodes. 

The last part of our study suggests that the derived general model can be easily applied 
to construct one- and two-dimensional lattices of active elements interacting by the nearest 
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neighbors type. Extension of the model to a three-dimensional case is straightforward. 

Finally, some limitations of our approach should be mentioned. First, the proposed 
model is not complete, there is no a time delay in pulse propagation among pacemakers, 
which can be very important for describing cardiac arrhythmias. Second, we represented 
cardiac tissue as a discrete one and used iterative approach to investigate its behavior. 
However, a large amount of realistic examples of active media is treated as continuous. 
Nevertheless, cardiac tissue is not a continuum, but is built up by discrete cardiomyocytes 
(or nodes with approximate dimensions 0.15 mm x 0.02 mm x 0.01 mm) [Kuramoto, 
1984]. 

Third and most important, to analyze the essential features governing dynamics of 
network of active elements, we have not included many important properties of the real 
conductive cells. These include the relaxation after stimulating, the prolonged (non-peak) 
form of pulses profiles, realistic topological structure, etc. Further investigations are re- 
quired to incorporate these features to the general combined model. 
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Figure 1: The model of two interacting pulse oscillators. 
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Figure 2: Model of mutually acting pulse oscillators. 
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Figure 3: Different types of tlie behavior of two bidirectionally interacting pacemakers. 
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Figure 4: Schematic representation of three pacemakers in the cardiac tissue. 
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Figure 5: Different types of the beliavior of tliree pulse oscillators. 
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Figure 6: 2D and ID lattices of coupled pulse oscillators. 
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